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Abstract. In this work we present a new simple but efficient scheme - Subsquares ap- 
proach - for development of algorithms for enclosing the solution set of overdetermined 
interval linear systems. We are going to show two algorithms based on this scheme 
and discuss their features. We start with a simple algorithm as a motivation, then we 
continue with a sequential algorithm. Both algorithms can be easily parallelized. The 
features of both algorithms will be discussed and numerically tested. 
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1 Introduction 

In this paper we address the problem of solving overdetermined interval linear 
systems (OILS). They can occur in many situations e.g. computing eigenvectors 
of interval matrices [2] or when solving various continuous CSP problems. There 
exist a lot of efficient methods for solving square interval linear systems. Solving 
overdetermined systems is a little bit more tricky, that is because we can not 
use some favourable properties of matrices like diagonal dominance, positive 
dennitcness etc. Nevertheless there are some methods - Rohn method [5], the 
least squares approach [3J, linear programming [3J, Gaussian elimination pQ or 
the method designed by Popova [7]. 

In our text [3J we showed that one of the best methods is using the least 
squares approach. This method returns a very narrow enclosure in a small time. 
But there is a problem, that the least squares always return solution, even if 
the system has none. That is sometimes not desirable. Other methods often 
rely on some kind of preconditioning which leads also to an overestimation and 
for some systems (e.g. with really wide intervals) can not be done. It is very 
difficult to develop one method suitable for every type of systems. We would 
like to present a scheme - Subsquares approach - which enables us to develop 
methods for solving overdetermined interval linear systems. Then we will move 
to a simple method and sequential method, both suitable for parallel computing. 
Before introducing the scheme and derived methods, it would be desirable to 
start with some basic interval notation and definitions first. 



2 Basics of interval arithmetics 



We will work here with real closed intervals x = [x, x] , where x < x. The 
numbers x, x are called the lower bound and upper bound respectively. 

We will use intervals as coefficients of matrices and vectors during our com- 
putations. The interval representation may be useful in many ways - it may 
represent uncertainty (e.g. lack of knowledge, damage of data), verification (e.g 
errors in measurement), computational errors (e.g. rounding errors in floating 
point arithmetic) etc. Intervals and interval vectors will be denoted in boldface 
i.e. x, b. Interval matrices will be denoted by bold capitals i.e. A, C . 

Another notion we will use is the midpoint of an interval x, it is defined 
as x c = (x + x)/2. By A c we will denote the midpoint matrix of A. When 
comparing two intervals we need the notion width of an interval x defined as 
w(x) = x — x. If u is an n-dimensional interval vector we will define "width" 
and "volume" of u as 

n n 

w(«) = w ( m *)' = n w ( u *)- 

i=l i=l 

The vector and matrix operations will be defined using the standard interval 
arithmetic, for definition of interval operations and further properties of the 
interval arithmetic do not hesitate to see e.g. [5]. 

We continue with definitions connected with interval linear systems. Let us 
have an interval linear system Ax = 6, where A is an m x n matrix. When m = n 
we will call it a square system. When m > n we will call it an overdetermined 
system. In the further text when talking about an overdetermined system, we 
will always use the notation Ax — 6, where A is m x n matrix. 

It is necessary to state what do we mean by the solution set of an interval 
linear system. It is the set 

£ = { x | Ax = b for some ie A, b <G b }. 

We shall point out, that this approach is different from the least squares method. 

If a system has no solution, we call it unsolvable. The interval hull is an 
n-dimensional box (aligned with axes) enclosing the solution set as tightly as 
possible. When we start using intervals in our computations (we have mentioned 
its advantage already), many problems become NP-hard [4]. So is the problem 
of finding the hull of the solution set [10]. It can be computed quite painfully 
using e.g. linear programming [3]. That is why we are looking for a little wider 
n-dimensional box containing the hull. The tighter the better. We call it interval 
enclosure. 

In this work we will provide numerical testing at various places in the text, 
therefore we rather mention its parameters here. The testing will be done on 
CPU Intel T2400, Core Duo, 1.83 GHz, with 2.50 GB memory. We used Matlab 
R2009a and toolbox for interval computation INTLAB v6 [TT] and Versoft vlO 
for verified interval linear programming. 
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All examples will be tested on random overdetermined systems. A random 
system is generated in the following way. First we generate random solvable 
point overdetermined system. Coefficients are taken uniformly from interval 
[—20, 20]. Then we inflate the coefficients of this systems to intervals of certain 
width. The original point system is not necessarily a midpoint system of the 
new interval system. Each of its coefficients is randomly shifted towards one of 
the bounds of an interval in which it lies. 

3 General Subsquares method scheme 

By a square subsystem (we will also call it a subsquare) of an overdetermined 
system we mean a system composed of some equations taken (without repeti- 
tion) from the original overdetermined system such that together they from a 
square system. Some of the possibilities are shown in the Figure [TJ For the 
sake of simplicity we will denote the square subsystem of Ax — b created by 
equations ii,i2, . . . ,i n as A{ il i2 = b{j x ,i 2 , ...,«„}• When we use some order 

(e.g. dictionary order) of systems (here it does not depend which one) the j-th 
system will be denoted Ajx = bj. 

Let us suppose we can solve a square interval system efficiently and quickly. 
We can for example take one of the following method - Jacobi method [5], 
Gauss-Seidel method [HIS], Krawczyk method etc. These methods usually 
can not be applied to overdetermined systems. Nevertheless, we can use the 
fact that we can solve the square systems efficiently together with the fact that 
the solution set of an overdetermined system must lie inside the solution set of 
its any subsquare. This follows from the fact that by adding new equations to 
the square system we can only make the solution set equal or smaller. 

When we chose some subsquares of an overdetermined system we can simply 
provide an intersection of their solution enclosures or provide some further work. 
We get a simple scheme for solving overdetermined interval linear systems - 
Subsquares Approach - consisting of two steps: 

1. Select certain amount of square subsystems of Ax = b 

2. Solve these subsystems and get together the enclosures 

If a method for solving OILS uses this kind of approach, we call it Subsquares 
method. As a motivation for this approach let us take the randomly generated 
interval system Ax = b (with rounded bounds), where 
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In the Figure[2]we can see the solution set and the hull of A^ 2 }X = b{i,2} ( re d 
color) and the same for A{ 2y 3}X — ^{2,3} (blue color). It can be seen that if 
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we provide intersection of the two hulls (or enclosures) the resulting enclosure 
might get remarkably tighter. 



a) b) c) 

Figure 1: Various square subsystems a), b), c) 




Figure 2: Enclosures and hulls of subsquares 
3.1 The simple algorithm 

If we solve square subsystems separately and then intersect resulting enclosures, 
we get a raw simple Algorithm [T] for solving OILS. 

Algorithm 1 Subsquares method - simple algorithm 
Require: A, b 

Ensure: enclosure x of the solution set of Ax = b 

x — [— oo, oo] n 

while not (terminal condition) do 

choose randomly a subsquare of Ax = b 
compute its enclosure x su ^ sq 

X . — X f~) X subsq 

end while 



This approach is a little bit naive, but it has its advantage. First, if we 
compute enclosures of all possible square subsystems, we get really close to the 
interval hull. The Table [T] shows the average ratios of widths and volumes of 
enclosures x su i, S(] , x ver returned by simple subsquares method and verifylss 
compared to the interval hull Xhull computed by linear programming. If we 
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have an m X n system, the number of all square subsystems is equal to (™) . 
However, we can see that for n small or for n close to m the number (™) might 
not be so large. That is why solving all subsquares pays off when systems are 
noodle-like or nearly-squared. 
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Table 1: Simple subsq. method solving all subsquares - enclosures comparison 

The second advantage is that Algorithm [l] can, in contrast to other methods, 
easily decide whether a system is unsolvable - if, in some iteration, the resulting 
enclosure is empty after intersection, then the overdetermined system is unsol- 
vable. The table [2] shows average number of random subsquares solved until the 
unsolvability is discovered (empty intersection occurs) . Each column represents 
systems of different coefficient radii. We can see that for systems with relatively 
small intervals unsolvability is revealed almost immediately. 
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Table 2: Simple subsq. method - unsolvability detection 



For most rectangular systems it is however not convenient to compute so- 
lutions of all or many subsystems. The choice of square subsystems and the 
solving algorithm can be modified to be more economical and efficient. 

3.2 The sequential algorithm 

When selecting subsquares randomly, they usually overlap. We can think of 
overlaps as a "meeting points" of square subsystems. They share some data 
(equations) there. That is why it would be advantageous to use this overlap 
to propagate a partial result of computation over one square subsystem into 
computations over other subsystems. When we are talking about immediate 
propagation of partially computed solution, our mind can easily come to Gauss- 
Seidel iterative method (GS). 
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This method starts with an initial enclosure a;' ). In fc-th iteration each entry 
of the current solution enclosure vector a;^ -1 ) might be narrowed using the 
formula 



1 ~ A- 



b t - (A a x[ k) + ... + A^^xfX + 



n x^\ 



Simply said, in each iteration this algorithm expresses x t from i-th equation. It 
uses newly computed values immediately. 

In our algorithm we will use GS iteration in a similar way for more square 
subsystems simultaneously. Again, we start with some initial enclosure x^ . In 
fc-th iteration we provide fc-th GS iteration step for all systems. The advantage of 
this technique is that we express each variable according to formulas originating 
from more systems. We expect the narrowing rate will be much better this way. 
Similarly as in simple GS, if in some point of computation empty intersection 
occurs, whole over determined system has no solution. 

Iterative methods usually require a preconditioning. We will use the pre- 
conditioning with A^ 1 . There are still two not answered problems yet - initial 
enclosure and terminal condition. To find x(°> , we can take the resulting en- 
closure of some other method. Or we can solve one square subsystem. The 
algorithm will terminate after fc-th iteration if e.g. 



(k) o-i) 
x\ - x\ 



< e and 



_(fe) _(fe- 
x ) —x) 



< e, 



for some small positive e and i = 1, . . . , n. 

Now we have to choose subsystems for sequential method. Before getting on, 
we shall consider the following desirable properties of a new algorithm inspired 
with four unfavorable features of the simple algorithm: 

1. We do not want to have too many square subsystems 

2. We want to cover the whole overdetermined system by subsystems 

3. The overlap of subsquares is not too low, not too high 

4. We take subsquares that narrow the resulting enclosure as much as possible 

First and second property can be solved by covering the system step by step 
using some overlap parameter. About third property, it proved itself to be a 
reasonable choice taking overlap « n/3. Property four is a difficult task to 
provide. We think deciding which systems to choose (in a favourable time) is 
still an area to be explored. Yet randomness will serve us well. 

Among many possibilities we tested, the following choice of subsystems 
worked well. During our algorithm we divide equations of overdetermined sys- 
tem into two sets - covered, which contains equations that are already contained 
in some subsystems, and waiting, which contains equations that are not covered 
yet. We also use a parameter overlap. 
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The first subsystem will be chosen randomly, other subsystems will be com- 
posed of overlap rows from covered and n — overlap rows from waiting. The 
last system is composed of all remaining uncovered rows and then some already 
covered rows are added to form a square system. This is described as Algorithm 
[2j The algorithm is not necessarily optimal, it should serve as an illustration. 
The procedure randsel(n, list) selects n random non-repeating numbers from 
list. The total number of subsquares chosen by this algorithm is 

m — n 
n — overlap 

The whole sequential algorithm is summarized as Algorithm [3j The function 
GS-iteration(Cx = d, y) applies one iteration of Gauss-Scidcl method on the 
subsquare Cx = d using y as initial enclosure . Method has-converged() returns 
true if terminal condition (e.g. the one mentioned earlier) is satisfied. 



Algorithm 2 Choosing square subsystems 
Require: A, b, overlap 
Ensure: set of subsquares of Ax = b 
systems 4— {set of square subsystems} 

covered 4— {numbers of covered equations by some subsystem} 
waiting 4— {1,2,..., m} {numbers of equations to be covered} 

while waiting 7^ do 
if covered = then 

indices 4— randsel(n, waiting) 

else if I waiting I < (n — overlap) then 

indices 4— waiting U randsel(n — jwaiting|, waiting) 

else 

indices 4— randsel(overlap, covered) U randsel(n-overlap, waiting) 
end if 

systems 4- systems U {Ai a dices2: = ^indices} 

covered <— covered U indices 

waiting 4— waiting \ indices 
end while 
return systems 



As we showed in [3], verifylss (using the least squares) from INTLAB is 
one of the best and quite general method for overdetermined systems that are 
solvable. That is why we wanted to compare our method with this method. 
During our testing we realized verifylss works fine with small intervals, how- 
ever it is not too efficient when the intervals become relatively large. We used 
enclosures returned by verifylss as inputs for subsquares sequential method 
and tested if our algorithm was able to narrow them. The Table [3] shows the 
results. Column rad shows radii of intervals of testing systems, we chose the 
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Algorithm 3 Subsquares method - sequential version 



Require: A,b,x m 

Ensure: enclosure x of the solution set Ax — b 
select square subsystems {Ait = bi, . . . , A^x — bk} 
x <- x<°) 

converged <— false 

while not converged do 

for i = 1 to k do 

x <— GS-iteration(Aia; = b iy x) 

end for 

converged «— has-converged() 
end while 
return x 



same radii for all coefficients of a system. We tested on 100 random systems. 
For each system we chose 100 random subsquares sets and applied the sequential 
algorithm on them. The fourth column shows average ratios of enclosure widths 
of x subsq an d x ver . Each random selection of subsquares set usually produces 
a different ratio of enclosure widths. For each system we chose one of the 100 
subsquares sets that produces the best ratio. The fifth column shows the aver- 
age value of the best ratios found for each of 100 random systems. Columns t ver 
and t su bsq show computation times (in seconds) of verif ylss and sequential 
method respectively. 
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Table 3: Subsquare method shaving the verifylss enclosure 



The larger interval radii are, the more intensively the sequential method 
sharpens verifylss enclosures. When an interval system has its midpoint 
system (A c x — b c ), or a system relatively close to it, solvable, we are able 
to compute the solution enclosure much tighter. This happens because of the 
preconditioning with A' 1 matrix. Selected results can be seen in the Table g 



8 



system 


overlap 


rad 


av. rat 


best av. rat 


15 x 10 


3 


0.1 


0.98 


0.89 


15 x 10 


3 


0.25 


0.85 


0.59 


15 x 10 


3 


0.35 


0.76 


0.40 


25 x 13 


5 


0.1 


0.99 


0.96 


25 x 13 


5 


0.25 


0.93 


0.74 


25 x 13 


5 


0.35 


0.76 


0.33 


37 x 20 


7 


0.1 


0.99 


0.97 


37 x 20 


7 


0.25 


0.89 


0.34 


50 x 35 


11 


0.1 


0.98 


0.82 



Table 4: Subsquares method shaving the verifylss enclosure (solvable midpoint 
system) 

3.3 Parallel algorithm 

The naive algorithm can be easily parallelized. All the square subsystems can 
be solved in parallel and then all enclosures are intersected. We have only one 
barrier at the end of computation. 

If we take a look at computation times in the Table [3j we realize verifylss 
is much faster. However, the time pay for gaining much precise enclosure using 
sequential method is not too high. Moreover, even the sequential algorithm can 
be parallelized. Propagation of newly computed enclosure can be guaranteed by 
sharing the enclosure vector x among processors as a global variable. If we use 
Jacobi formula instead of Gauss-Scidcl formula for one iteration, the computa- 
tion becomes of kind SIMD - single instruction multiple data. Therefore it could 
be used even on GPUs - one pipeline narrows one variable from the interval en- 
closure vector, a bunch of pipelines computes over one subsquare. Nevertheless, 
shared vector x might create a bottleneck. We believe this could by prevented 
by the following behaviour of each pipeline. When reading, each pipeline does 
not lock corresponding shared variable. After each iteration it overwrites a 
shared variable only if it has improved the value currently stored there. This is 
inspired with an observation, that in one iteration not many computations over 
different subsquares improve the same variable. However, we assume that there 
exist much more efficient ways how to make parallel subsquares methods more 
efficient and memory collision avoiding. 

4 Conclusion 

In this paper we introduced a simple but efficient scheme - subsquares approach 
- for enclosing the solution set of overdetermined interval linear systems. The 
first method derived from this scheme was a little bit naive, but for noddle-like 
or nearly-square systems it was able to find almost the interval hull. The second 
method was a little bit more sophisticated but still quite simple. It worked well 
on interval systems which coefficients were composed of wide intervals. This 
method was able to significantly sharpen enclosures produced by verifylss. 
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Both methods were able to detect unsolvability of OILS. Moreover they could 
be easily parallelized. In the second method we chose the square subsystems 
randomly, that is why sharpening produced by this method had variable results. 
There is a question whether for each OILS there exists a deterministically cho- 
sen set of subsquares which gives the best possible enclosure, so we can avoid 
the randomization. 
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